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THE  TRANS  IATIONAL  VELOCITY  OF  SURFACE  SHIPS  AND  SUBMARINES:  A  CCMPOTHt  PROGRAM 

This  report  is  part  of  a  continuing  study  of  the  Interaction  of  the  underwater 
explosion  shock  rare  with  the  ocean  bat  ton.  The  computer  program  described  In 
this  paper  was  primarily  written  to  calculate  the  translational  Telocity  Induced 
in  surface  ships  by  bottom  reflected  shock  waves.  These  calculations  provide  a 
method  of  comparing  the  damage  producing  potential  of  the  reflections  for  various 
bottom  materials.  The  work  was  done  under  the  supervision  and  cooperation  of 
Dr.  n.  G.  Snay  (240). 

This  study  was  supported  by  the  Defense  Atomic  Support  Agency  through  Task 
DASA-HA  002-20  P.  106  (Energy  Focussing  and  Refraction  Effects). 
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THE  PEAK  TRANS  IATIONAL  VELOCITY  OF  SURFACE 
SHIPS  ANT  SU94AKNES:  A  COMPUTER  IROQRAM 

1.  INTRODUCTION 

The  peak  translational  velocity  (PTV)  of  the  center  of  gravity  of  a  naval 
ship  or  submarine  induced  by  underwater  explosion  shock  vaves  is  generally  used 
to  describe  the  degree  of  impairment  of  their  mobility  and  weapon  delivery  capa¬ 
bilities.  The  model  presently  being  used  to  calculate  the  PTV  is  that  developed 
primarily  for  submarines  by  W.  W.  Murray  (reference  (l)).  This  model  treats  the 
interaction  of  an  exponentially  decaying  acoustic  plane  wave  with  an  infinitely 
long  cylinder.  For  pulses  of  nuclear  dimensions  the  assumption  of  plane  inci¬ 
dent  waves  is  usually  Justified  because  the  ranges  considered  are  large  compared  to 
the  dimensions  of  the  ship  or  submarine. 

In  the  application  of  Murray's  theory  to  waves  vtoich  have  been  reflected 
from  the  ocean  bottom  or  refracted  by  velocity  gradients  in  the  ocean  one  encounters 
the  need  for  calculating  the  PTV  for  wave  shapes  other  than  exponential.  One  of 
the  best  ways  to  make  such  a  calculation  for  an  arbitrary  wave  shape  is  through  a 
superposition  of  step  wave  responses.  Murray  has  calculated  the  step  wave  trans¬ 
lational  velocity  curve  and  also  the  step  wave  acceleration.  J.  A.  Goertner 
(in  a  confidential  report)  has  written  a  computer  program  which  uses  Murray's 
curves  to  calculate  the  PTV  for  an  arbitrary  incident  wave  by  decomposing  the 
wave  into  a  sum  of  step  waves.  This  program  has  been  used  successfully  in  cal¬ 
culating  the  PTV  of  refracted  waves,  but  is  not  well  suited  for  bottom  reflection 
studies . 

In  this  paper  Murray's  theory  is  described  briefly,  and  a  computer  program  Is 
explained  -Which  computes  the  PTV  for  an  arbitrary  wave  shape  in  a  somewhat  dif¬ 
ferent  manner  than  Goertner' s  program.  The  incident  pulse  used  in  the  program  of 
this  paper  may  have  a  singularity  of  the  logarithmic  type  such  as  encountered  in 
supercritical  bottom  reflections.  The  PTV  is  calculated  by  a  convolutioi  integral 
containing  the  incident  wave  shape  and  the  step  wave  acceleration.  The  curve  of 
the  step  wave  acceleration  has  been  recalculated  so  that  the  model  can  be  more 
closely  followed  than  is  possible  using  Murray's  curve.  The  theory  is  extended 
to  surface  ships,  and  the  program  calculates  the  PTV  for  both  surfaced  and  sub¬ 
merged  targets. 
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2.  THEORY  FOR  CALCULATING  THE  TRANS IATIONAL 
VELOCITY  OF  AN  INFINITELY  LONG  CYLINDER 


2.1  Assumptions 

Murray  derived  his  equations  for  a  rigid  and  neutrally  buoyant  cylinder  of 
radius  a.  It  Is  assumed  that  the  displacement  of  the  cylinder  from  its  initial 
position  is  small  compared  to  its  radius.  The  equations  vere  derived  for  athwart- 
ship  attack;  that  is,  the  wave  front  is  parallel  to  the  longitudinal  axis  of  the 
cylinder. 

2.2  Translational  Velocity  of  a  Submerged  Cylinder 
Let  the  incident  vave  be  given  by 

p(t)  *  Pp  exp  [-  (t  -  R/cv)/G]  ,  t  >  R/cy  (2.2.1) 

where  t  is  the  time,  R,  the  distance  from  the  source  to  the  target,  G,  the  time 
constant  of  the  exponential  shock  wave  (usually  denoted  by  0),  and  c^,  the  sound 
velocity  of  water.  The  peak  pressure  of  the  wave  is  p^.  For  this  exponential 
pulse  Murray  obtained  the  following  equation  for  the  translational  velocity  of  a 
totally  submerged  cylinder 


v(i) 


i«o+v 

/ 

-i'T'+V 


[*lT.-U3_  at  , 

z  (z+q)  Kg  (z) 


(2.2.2) 


where  the  integration  variable  z  is  a  complex  magnitude  and  0^  is  the  density  of 
water.  Hie  symbol  T  denotes  the  reduced  time  T  »  c^  t/a,  and  q  is  the  reduced 
radius  q  *  a/cjG.  For  •>  step  wave  G  becomes  infinite,  and  we  have  q  =»  0.  The 
function  Kg(z)  is  the  modified  Bessel  Function  of  the  second  kind  of  the  order 
two.  The  path  of  integration  is  to  be  taken  in  the  right  half  of  the  complex 
plane,  hence  the  constant  v  must  be  real  and  positive.  For  practical  purposes,  a 
good  choice  of  v  is  unity. 

2.3  Reduced  Gtep  Wave  Acceleration 

Upon  differentiating  v(t)  and  setting  q  ■  0,  the  desired  expression  for  the 
reduced  step  wave  acceleration  of  the  cylinder  is 
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The  reduced  step  wave  acceleration  A(t)  plays  the  role  of  a  Green's  function 
for  the  problem.  The  translational  velocity  V(t)  from  an  arbitrary  incident  wave 
P(t)  can  be  written 

v(T)  “  7T"  /  p(Wc  )  a(t  -  q)  dq  ,  (2.4.1) 

pv  w 

o 

where  T  =  c^t/a.  If  the  integration  variable  is  changed  so  that  it  has  the 
dimensions  of  time,  V(t)  is  then  given  by 

t 

V(t)  *  f  p(u)  a(t  -  c  u/a)  du  .  (2.4.2) 

°wa  J  w 

o 
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This  1b  the  equation  used  to  calculate  V(t)  in  the  FTV  PROGRAM  described  in 
Section  3» 

2.5  Translational  Velocity  of  a  Surface  Ship 

To  apply  the  above  equation  to  the  response  of  a  surface  ship,  two  assumptions 
are  made:  (l)  the  target  is  considered  to  be  a  cylinder  floating  on  the  surface 
with  its  axis  at  the  water  line.  (2)  the  vertical  translational  velocity  is  as¬ 
sumed  to  be  twice  the  vertical  component  of  the  translational  velocity  the  cylin¬ 
der  would  acquire  deeply  submerged.  The  horizontal  motion  of  the  ship  is  not 
taken  into  account. 

These  assumptions  are  usually  made  for  calculation  of  damage  to  surface 
ships,  although  it  is  realized  that  it  may  be  an  oversimplification.  Effects 
such  as  cavitation  are  also  ignored.  This  process  is  known  to  occur  below  ships 
and  may  be  of  importance. 

Under  the  above  assumptions,  the  vertical  translational  velocity  of  a  float¬ 
ing  cylinder  when  subjected  to  a  pressure  pulse  p(t)  is  then 

Vs(t)  -  2V(t)  cos  n.  ,  (2.5.1) 

and  t 

vs(t)  =  f-f22-®  /  P(u)  A(t  -  cvu/a)  du  ,  (2.5.2) 

w  o 

where  a  is  the  angle  between  the  plane  wave  front  and  a  normal  to  the  water  sur¬ 
face,  sometimes  called  the  incident  angle. 

Hie  program  described  in  Section  3  calculates  both  V(t)  and  V  (t)  and  their 

8 

maximum  values,  tk^  peak  translational  velocities  FTV. 

2.6  Comparison  of  the  Responses  of  a  Target  to  Exponential  an  Supercritical 
Bottom  Reflected  Pulses. 

The  experimental  data  correlating  the  shock  damage  from  an  unde,  water  explo¬ 
sion  to  the  peak  translational  velocity,  PTV,  have  been  obtained  for  free  water 
pulses  or  for  free  water  pulses  cut  off  by  surface  reflections.  Both  of  these 
pulse  shapes  are  initially  exponential.  Pulse  shapes  encountered  in  the  study  of 
supercritical  bottom  reflections  are  not  exponentials,  and  the  question  arises 
whether  the  same  shock  damage  results  if  the  FTV's  sire  the  same.  Two  examples  of 
these  ptlses,  along  with  an  exponential,  are  given  in  Figure  2.6.1.  As  shown  in 
Figure  2.6.2  these  pulses  produce  the  same  FTV  on  a  cylinder  of  radius  22  ft. 
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FIG.  2.6.1  FREE  WATER  AND  BOTTOM  REFLECTED  PULSES  PRODUCED  BY  A  10  KT  CHARGE 


TRANSLATIONAL  VELOCITY  V/PTV 
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TIME  (MSEC) 


FIG.  2.6.2  RESPONSES  OF  A  CYLINDER  OF  RADIUS  22  FT.  TO  FREE  WATER  AND 
BOTTOM  REFLECTED  PULSES  PRODUCED  BY  A  10  KT  CHARGE 
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Ignoring  the  early  parts  of  the  responses  to  the  reflected  pulses,  these  curves 
have  roughly  the  same  shape  around  the  peak  as  the  response  to  the  exponential. 
After  the  peak  the  cut-off  exponential  response  deviates  much  more  than  those  of 
the  bottom  reflections.  Having  the  same  PTV  and  similar  accelerations,  the  pulses 
of  Figure  2.6.1  are  expected  to  cause  the  same  degree  of  damage.  This  means  that 
PTV  damage  criteria  derived  for  exponential  pulses  can  also  be  applied  to  super¬ 
critical  bottom  reflections  and  other  similar  non-exponential  pulses. 


3.  COMPUTER  PROGRAM  FOR  CALCULATING 
PEAK  TRANS IATIONAL  VELOCITY 


3*1  General  Program  Description 

The  peak  translational  velocity  program  or  simply  the  PTV  PROGRAM  has  been 
written  in  FORTRAN  IV  for  the  NOL  CJC  6400  computer.  A  complete  listing  is  given 
in  Appendix  C.  This  program  calculates  the  PTV  for  both  surface  ships  and  sub¬ 
marines  using  the  theory  described  in  Section  2. 

The  PTV  PROGRAM  is  composed  of  seven  subroutines:  PTV,  FV,  FI,  XMAX,  VTAB, 
PTAB,  and  FGI.  The  package  is  us  d  by  calling  subroutine  PTV  from  a  main  or 
executive  program  written  by  the  user  vfaich  supplies  the  pressure  time  history  p(t). 
The  PTV  is  obtained  from  equations  (2.4.2)  and  2. 5*1).  But  in  order  that  we 

may  integrate  numerically  over  a  singularity  in  p(t)  at  t  *  t  ,  or  t  =  t  =  c  t  /a, 

Ic  c  w  c 

t-t  ,the  integration  variable  u  is  changed  as  follows: 
c 


for  u  i  t 

w  «  (t  -  u)1//2 

c 

c 

for  u  ?  t 

c 

z  =  (u  -  t  )1/2 
c 

Equation  (2.4.2)  then  becomes  for  t  >  t 

c 


o  (  v(t  ) 

V(t)=  -  <  -  f c  p(u)  A(t  -  c  u/a)w  dw 

°wa  ‘  1  x  w 

W  w(0) 

*(t ) 

+  f  p(u)  A(t  -  cyu/a)z  dz  l  (3.1.1) 

z(tc) 
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where  t  »  c^t/a.  These  integrals  are  evaluated  in  FUNCTION  FV  using  the  Gaussian 

quadrature  of  FUNCTION  FGI.  From  V(t)  we  then  obtain  V  (t)  using  equation  (2.5*1). 

■ 

The  A(t)  curve  which  has  been  calculated  by  the  procedure  of  Appendix  A  is 

stored  in  the  arrays  QQX  and  QQY  in  FUNCTION  FI.  lhe  reduced  time  t  is  in  QQX  and 

A  in  QQY.  The  function  A(t  -  c  u/a)  is  evaluated  from  these  arrays  by  quadratic 

interpolation  in  FUNCTION  VTAB.  Similarly  p(u)  is  determined  by  interpolation  in 

VTAB  of  the  arrays  QX  and  QY  which  hold  the  time  t  in  seconds  and  the  incident 

pressure  in  psi.  Near  the  singularity  at  t  *  t  the  FUNCTION  PEAB  performs  the 

c 

quadratic  interpolation  for  the  pressure. 

The  convergence  of  the  Integrals  in  equation  (3.1.1)  is  made  possible  because 

lim  w  ln|t  -  t  |  »  lim  z  ln|t-t  I  =  0.  (3.1.2) 

u-*t  c  u  -  t  lcl 

c  c 

As  Implied  in  equation  (3.1.1)  the  variables  v  and  z  are  used  for  integration 

over  the  whole  range  of  t.  Little  difficulty  is  encountered  in  the  numerical 

integration  if  the  pressure  pulse  p(t)  has  no  rapidly  changing,  high  amplitude 

contributions  far  from  the  peak  at  t  *  c  t  /a. 

C  V  C 

The  values  of  V(t)  and  V  (t)  depend  on  the  previous  pressure  history.  Since 

8 

A(t  )  is  very  small  for  t  >  8,  the  integration  range  is  restricted  to  at  most  from 
u  *  t  -  8  to  u  *  t,  Thus  if  significant  rapidly  changing  pulses  occur  away  from 

t  by  about  t  «  8,  the  FTV  PROGRAM  can  be  applied  to  each  peak  separately  since 

c 

the  target  response  from  one  pulse  is  essentially  damped  out  before  the  arrival  of 

the  next  pulse.  The  actual  PTV  can  then  be  found  from  the  maximum  of  these  results. 

The  maximum  or  peak  values  of  V(t)  and  V  (t),  the  PTV' s» are  obtained  as  fol- 

8 

lows.  An  initial  search  for  a  maximum  velocity  is  made  from  some  t  *  t  to  t  «  t^. 
The  values  of  tQ,  t^,  and  the  number  of  steps  are  prescribed  by  the  user  in  the 
call  to  subroutine  FTV.  Then  several  Iterations  are  made  around  this  maximum. 
Subroutine  XMAX  determines  the  maximum  value  of  the  translational  velocity,  but 
subroutine  FTV  controls  the  iteration  and  makes  the  calls  to  FUNCTION  FV  i&lch 
sets  up  the  integration  for  V(t).  Iteration  terminates  irtien  the  relative  dif¬ 
ference  between  the  two  largest  absolute  values  of  V(t)  is  less  than  .001.  If  the 
iteration  does  not  converge  after  five  cycles,  iteration  is  also  terminated  and  a 
warning  is  printed.  In  either  case  values  of  the  PTV  for  submerged  targets,  the 

maximum  of  V(t),  and  for  floating  targets,  the  maximum  of  V  (t),  are  returned  to 

8 

the  main  program. 
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3.2  Use  of  the  PTV  ffiOGRAM 

To  use  the  PTV  PROGRAM  subroutine  package  a  main  program  must  be  set  up  by 
the  user  to  supply  the  incident  pulse  p(t).  The  time  in  seconds  and  the  pressure 
in  psi  must  be  stored  in  the  arrays  QX  and  QY  as  mentioned  previously  in  Section 
3*1*  When  the  pressure  history  is  short  compared  to  the  target  transit  time  a/ cw, 
the  PTV  is  likely  to  occur  at  a  time  beyond  the  last  value  of  the  pressure  history. 
Thus  to  provide  for  extrapolation  beyond  the  end  of  the  actual  pressure  history 
the  first  unused  storage  of  the  QX  array  should  be  set  to  some  very  large  value 
as  1.GE20.  The  corresponding  QY  storage  should  be  set  to  zero  or  some  other 
appropriate  asympto'  ic  value  of  p(t). 

The  QX  and  QY  arrays  are  transferred  to  the  PTV  PROGRAM  by  COMMON  storage. 

The  statements  COMMON  /QXY/QX,QY  and  DIMENSION  QX(lOOO),  QY(lOOO)  must  be  in  the 
main  program.  In  t ubroutines  PTV  and  71  the  additional  coranon  storage  is  used: 
COMMON/QIs/lS .  This  statement  is  not  needed  in  the  main  program. 

Once  the  pressure  history  has  been  defined,  the  peak  translational  velocity 
is  then  obtained  by  calling  subroutine  PTV  as  follows: 

CALL  PTV  (TIMER2,  T3,  T4,  T5,  HAD,  PTS,  OPTION,  COSA,  RHOW,  CWAT,  T,  V,  VS). 

INPUT  The  following  variables  are  inputs  to  subroutine  PTV: 

TIMER2  Time  t  in  seconds  of  the  singularity  or  peak  of  the  incident 

c 

pulse.  For  a  simple  exponential  pulse  set  TIMER2  =  0.  The  pres¬ 
sure  at  a  singularity  should  be  set  to  some  number  with  absolute 
value  greater  than  1.0E20  as  a  signal  to  the  interpolation 
subroutine  FTAB. 

T3  Signals  the  approach  of  the  singularity  of  the  incident  pulse 

p(t).  If  there  is  no  singularity  set  T3  =  TIMER2.  When  there 
is  a  singularity,  T3  should  Have  a  value  such  that  there  are  in¬ 
cluded  at  least  two  points  of  the  QX  array  on  each  side  of  the 

singularity  in  the  time  interval  T3  ^  t  <  2t  -  T3. 

c 

T4  Smallest  time  t  at  which  the  translational  velocity  is  to  be 

calculated.  If  the  peak  of  p(t)  occurs  at  or  near  zero,  use 
T4  «  0.  In  other  cases  T4  (and  T5  below)  can  be  determined  by 
remembering  that  the  translational  velocity  at  time  t  is  calcu¬ 
lated  using  the  pressure  history  of  the  interval  t  -  8a/cy  to  t. 
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T5  Largest  time  t^  at  which  the  translational  velocity  is  to  be 

calculated. 

RAD  The  cylinder  radius  a  in  feet. 

PTS  The  number  of  times  at  which  the  translational  velocity  is  to  be 

calculated  in  the  initial  search  for  the  HTV.  This  search  is 
made  in  the  time  interval  t  <  T5-  The  maximum  value  PTo 
can  be  is  50. 

OPTION  Controls  printing  in  subroutine  PTV.  There  is  no  printing  if 

OPTION  >  0.  There  is  printing  if  OPTION  »:  0. 

COSA  cos  a.  See  Section  2.5  for  an  explanation  of  a. 

RHOV  Density  of  water  in  gm/em^. 

CWAT  Sound  velocity  of  water  c^  in  ft/sec. 

OUTPUT  T'  *  following  variables  are  outputs  returned  to  the  main  program.  When 
OPTIC  ,  these  results  are  printed  out  in  subroutine  PTV. 

T  Time  t  in  seconds.  The  time  of  the  PTV  is  returned  to  the  main 

program. 

V  The  translational  velocity  V(t)  in  ft/sec  of  a  submerged  target. 

The  PTV  is  returned. 

VS  The  translational  velocity  V  (t)  in  ft/sec  of  a  floating  *arget. 

s 

The  PTV  is  returned. 

A  sample  print  out  for  a  pressure  pulse  p(t)  =  exp  (  -  125t),  or  p(t)  = 
exp  (-  t)  when  a/c^  =  .008,  is  shown  in  Appendix  D.  The  input  to  subroutine  PTV 
is  included  in  the  print  out. 

3.3  Important  FORTRAN  Symbols  Not  Included  in  the  Call  to  Subroutine  PTV 
Dimensioned  Variables 
SUBROUTINE  PTV 

QX,  QY  Time  in  seconds  and  pressure  in  psi  of  the  incident  pulse. 

These  arrays  must  be  defined  in  the  user's  executive  program. 

QQX,  QQY  Reduced  time  and  reduced  acceleration  of  step  wave.  These  arrays 

appear  in  FUNCTION  FI. 
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IS(1) 

IS(2) 

G 

G(l)  =  T 
G(2) 

G(3) 

G(M 

G(5) 

0(6) 

A,  C 

Non-  jensloned 
SUBROUTINE  PTV 
UT 

M 

Tl,  T2,  VI,  V2 
VS1 

FUNCTION  FV 
N  -  18 


NN,  Nl,  NNN 


X 


Index  for  the  beginning  of  the  interpolation  search  in  QQX  array. 

Index  for  the  beginning  of  the  interpolation  search  in  QX  array. 

Array  for  transferring  to  FUNCTIONS  FV  and  FI  variables  in  the 
integrand  of  V. 

Time  t 


t 

c 


Signal  for  FUNCTION  FI.  In  equation  (3.1.1)  G(4)  *  -  1.0  for 
the  first  integral  and  +  1.0  for  the  second  integral. 


<te  -  T3) 


v(0) 


1/2 


t1/2 


Storage  for  time  and  V(t).  Used  by  subroutine  XMAX  to  determine 


the  maximum 


V(t) 


C(M)  and  the  next  largest  value  C(Ml). 


Variables 


Incremsnt  of  time. 

See  A  and  C  above. 

Temporary  storages  of  A(M),  A(Ml),  C(m),  c(Ml). 

Value  of  V  (t)  vhen  V(t)  -  C(Ml). 

8 

Tbe  integrations  of  equation  (3.1.1)  are  performed  using  a  four 
point  Gaussian  quadrature  per  subinterval  of  integration.  N  is 
the  maximum  number  of  subintervals  allowed  for  the  total  inte¬ 
gration  Interval. 

The  number  of  subintervals  of  integration  used.  NN  is  used  if 

the  total  integration  Interval  does  not  Include  t  .  Nl  and  NNN 

c 

are  used  if  t  is  included:  Nl  for  the  Integration  variable 
c 

u  t  and  NNN  for  u  >  t  . 
c  c 

t  -  8c ^/a  used  to  restrict  integration  to  the  interval 
t  -  8  to  t. 
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Zl,  Z 2,  Z3 


FV 

FUNCTION  FI 

Z 

X 

XD 

P 

FI 


Limits  of  Integration  in  equation  (3.1.1).  In  the  calls  to 
FUNCTION  FGI  the  first  variable  is  the  lower  limit  of  Integra 
tion,  the  second  is  the  upper  limit. 

The  sum  of  the  integrals  of  equation  (3*l.l). 

Integration  variables  v  and  z. 

Time  corresponding  to  integration  variable  u. 

Reduced  time  equal  to  c  •  (t  -  u)/a. 

Interpolated  pressure  at  time  X. 

Integrands  of  equation  (3.1.1). 
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APPENDIX  A 

CALCULATION  OP  THE  REDUCED  STEP  WAVE  ACCELERATION  A(t) 


In  order  to  evaluate  A(t)  from  equation  (2.3*1)  it  is  necessary  to  transform 
the  integral  in  the  complex  plane  to  a  real  integral.  Murray  has  accomplished 
this  transformation  by  using  a  series  expansion  ihe n  t  is  small,  up  to  about  t  *  1, 
and  contour  integration  at  larger  values  of  t.  However,  to  obtain  a  more  accurate 
A(t),  ve  have  used  the  more  direct  approach  explained  below. 

The  integration  variable  z  can  be  written  z  *  x  +  iy  for  x  and  y  real.  If 
the  integration  path  is  taken  along  the  line  x  *  y  *=  1,  z  becomes  z  =  1  +  iy.  The 
complex  functions  in  the  integrand  of  A(t)  can  then  be  separated  into  their  real 
and  imaginary  parts: 

t2  »  (l  -  y2)  +  i2y  , 

exp  [z(t-1)]  *  exp  (t-1)  cos  [y(T-l)]  +  i  exp  (t-1)  sin  [y(T-l)]  , 

and  K2(z)  =  Re(K2)  +  i  Im(K2)  , 

where  Re(K2)  and  Im(K2)  denote  the  real  pert  and  the  imaginary  part  of  K2(z). 
Explicit  expressions  from  which  Re(Kg)  and  Im(K2)  can  be  obtained  will  be  given 
later.  On  substituting  the  above  functions  in  A(t),  equation  (2.3»l),  and  then 
separating  real  and  imaginary  parts  of  the  integrals  one  obtains 


A(l)  .  «*>  /  h  coe  *  *2  8ln  [y(T'1)]  ^ 

^  '  19  X  W  ^ 


E1  E2 


where 


+  i  f  Ei 8ln  [y(T-i)]  -  e2  oo*  [y(T-D]  ^  , 


T  2 

E1  +E2 


(A.l) 


\  -  (1  -  y2)  Re(K2)  -  2y  Im(K2) 


(A.2) 


and 

Eg  *  (l  -  y2)  Im(K2)  +  2y  Re(K2)  .  (A. 3) 

A  substitution  of  -y  for  y  in  equation  (A.l)  shows  that  the  integrand  of  the 
first  integral  of  A(t)  is  even  and  the  integrand  of  the  second  integral  is  odd. 


A-l 
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Hence  the  second  Integral  is  zero  and  A(t)  is  a  real  function  vhich  can  be 
written 


A(t)  -  *  exp  (x-1)  /  bi  •  •  (A.4) 

'  '  ^  ^  j.  n  ^ 


V  +  E2 


For  y  <  15  we  have  calculated  Kp(z)  from  the  expression 

GO 

Kg(z)  =  J*  exp  (  -  z  cosh  O  cosh  2$  d§  . 


(A.  5) 


Separating  the  exponential  into  its  real  and  imaginary  parts  and  substituting 

z  ■  1  +  iy,  ve  obtain 
00 

K^(z)  =  I  exp  (  -  cosh  $)  cos  (y  cosh  $)  cosh  2i  d $ 
o 

GO 

-  i  j  exp  (  -  cosh  t )  sin  (y  cosh  i)  cosh  2$  d$  .  (A. 6) 

o 

Substitution  of  this  expression  for  Kg(z)  into  equations  (A. 2)  and  (A. 3)  yields 
the  following  expressions  for  and  Eg 


CO 


El- 

( 

e 

1*  _  1 
'  (1-y2)  u  +  2y  zj  ux  d$ 

> 

0 

(A.7) 

E2  "  J 

C 

1  |_2y  u  -  (l-y2)zj  Ux  d* 

> 

(A.8) 

*Aiere:  U  =  cos  (y  cosh  $) 

Z  =  sin  (y  cosh  i) 

*  exp  (-  cosh  t)  cosh  2$. 

These  integrals  converge  very  rapidly  because  of  the  expression  U^which  approaches 

zero  like  exp  (-  exp  $)  for  about  $  =  4  or  larger.  It  is  shown  in  Appendix  B  that 

the  error  in  truncating  the  integration  in  E.  and  E  at  $  =4.5  is  less  than  1 
13  12 

part  in  10  . 
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Even  though  the  integrals  converge  rapidly,  they  become  increasing  more 
difficult  to  evaluate  numerically  ns  y  Increases  because  of  the  oscillatory  fac¬ 
tors  U  and  Z.  At  about  y  *  15  an  asymptotic  expansion  for  evaluating  K„(z) 
becomes  more  practical. 

For  y  between  15  and  1000,  the  following  asymptotic  expansion  (reference  (?)) 
of  K2(z)  is  used 


*,<■>  -  <s>1/2  -  (->  IviHfc*  *  -I  • 


21  (8z Y 


(A. 9) 


where  again  z  ■  1  +  iy.  Near  y  *  15  nine  terms  of  the  series  in  brackets  are 

used,  i.e.,  the  lowest  ordered  terra  used  is  of  the  order  l/z  .  Retaining  nine 

-9 

terms  Insures  that  the  series  truncation  error  for  y  *  15  is  less  than  2  x  10  . 

Between  y  ■  15  and  y  =  1000  fewer  terms  are  needed  for  larger  y;  however,  a  suf¬ 
ficient  number  of  terms  are  retained  so  that  the  truncation  error  is  less  than 
that  at  y  ■  15 . 

The  Integral  for  A(t)  from  y  *  1000  to  infinity  is  calculated  from  an  approxi¬ 
mate  equation  obtained  by  neglecting  terms  of  order  l/y2  or  smaller  compered  to 
one.  From  equation  (A. 9)  the  approximate  relation  for  K^d  +  iy)  is  obtained 

l/g 

Kg( l+iy )  *  UJ  cxp  [C06  ('y+*'  -  1  sin  (y+t)]  [1  -  (A. 10) 


where 


cos 


* 


(l  +  ~)  and  sin  ♦ 

2y 


(A.n) 


Substituting  the  real  and  imaginary  parts  of  K  from  equation  (A.IO)  into  equations 
(A. 2)  and  (A. 3)  and  neglecting  terms  of  order  l/y  ,  the  following  relations  are 


obtained : 


Ea  = 


l/2 

i/a 


exp  (-1)  [-yL  cos  (y  +  t>)  +  y  sin  (y  +  *)] 
exp  (-1)  [y2  sin  (y  +  *)  +  ^  y  con  (y  +  ifc)] 


f 


(A. IP) 
(A. 13) 


and 

Ex2  +  E22  y3  exp  (-2)  . 
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Combining  the  above  equations  ulth  equation  (A. 4),  the  remainder  R(a)  of  the  A(t) 
integral  from  y  =  a  to  infinity  is 


r  2V*/2 

R(a)  «  -I  exp  (t) 


r”3/i? 


^  3 in  (Ty  +  ii)  -  cos  (ry  +  +  )J  dy  ,  (A.15) 


where  the  trlgnometric  relations  for  the  sine  and  cosine  of  the  sum  of  two  angles 
have  been  used  and  where  the  lower  ordered  terms  have  been  neglected.  For  the 
numerical  computations  a  *  1000  is  used. 

Using  similar  manipulations  as  above  and  substituting  cos  $  and  sin  t  from 
equation  (A.ll),  R(a)  can  be  written 

R(a)  exp  (t )  j  y“3^2  [(l  +  |£)  sin  Ty  -  (l  -  |^)  cos  Ty]  dy.  (A.l6) 

a 

Integration  by  parts  can  then  be  used  to  obtain 

R(»)  **  “  9T)  [a  (sin  aT  "  cos  aT)  +  ^ 2ttt  (l  “  s(aT) 

-  C(aT ) )]  +  2  e~3/2  £gin  aT  +  cos  aT]  |  f  (A. 17) 


\toere  S(aT)  and  c(aT )  are  commonly  called  Fresnel's  integrals  and  are  defined 

aT  aT 


s(aT )  *  —  r  SiS-JS  dx  and  c(aT)  ,  _1_  |  cos_x  dx  # 

JT  m  .  JT 


(a. 18) 


These  integrals  have  the  following  asymptotic  expansions  (see  reference  (3)): 


c(z)  «  \  +  fi 
2  /2rrz  L 


(2z  f  (2«r 


•  •  • 


] 


cos  z 


i&l 

(SO3 


+ 


•  •  • 


s(z) 


cos  z 


1 


-111.  ♦  .vjl-jl-j 

Uzf  (Lr 


(A.  19) 


A-l* 


J  2rrz 


NOLTR  71-65 


In  summary,  to  evaluate  A(t)  from  equation  (A.  4)  and  E^,  which  are  defined 
by  equations  (A. 2)  and  (A. 3)#  are  given  by  equations  (A. 7)  and  (A. 8)  for  0  <  y  «  15 
and  obtained  from  equation  (A. 9)  far  15  s  y  s  1000.  The  Integral  from  y  =  1000  to 
Infinity  Is  R(a),  equation  (A. 21),  where  a  3  1000.  The  A(t)  curve  shown  In 
figure  2.3*1  was  calculated  by  the  above  method  on  the  NOL  191  7090  computer.  A 
table  giving  the  A(t)  array  to  six  decimal  places  is  contained  In  the  DATA  state¬ 
ment  of  FUNCTION  FI  of  the  FORTRAN  listing  of  the  PTV  PROGRAM  which  Is  given  In 
Appendix  C. 
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APPENDIX  B 

CONVERGENCE  OF  THE  INTEGRALS  AND  E£ 

It  is  the  object  of  this  section  to  shov  that  the  improper  integrals,  equa¬ 
tions  (A. 7)  and  (A.8j,  used  to  obtain  and  are  convergent.  We  also  obtain 
an  upper  bound  on  the  error  introduced  by  stopping  the  integration  to  infinity  at 
a  finite  value  of  the  Integration  variable 

It  is  commonly  proved  in  text  books  of  integral  calculus  that  improper  inte¬ 
grals  of  the  form  of  E1  and  E  are  convergent  if  the  integral  of  the  absolute 
value  of  the  integrand  is  convergent.  The  converse  does  not  necessarily  hold. 
Denote  the  integrand  of  E^  by  E^'  and  that  of  E^  by  E^ ' .  Since  in  general 
|a±b|^|a|+|b|,  j  sin  aj  >  1,  |  cos  a  |  si,  and  |  ab  j  =  |a|  |b| ;  ve  obtain 

|  y|  -  |[(1  -  y2)  U  ♦  2y  2]  ux  | 

<  (l  +  y2  +  2y)  exp  (-  cosh  $)  cosh  2$.  (B.l) 

This  result  also  holds  for  E^'  .  Since 

cosh  2$  ■  [exp  (2$)  +  exp  (-2*)]/2  "  exp  (2$), 
the  above  inequality  can  be  simplified  to 

|  E^  j  <  (l  +  y2  +  2y)  exp  (2$  -  cosh  l)  . 

At  ♦  *  4.5,  cosh  $  **  l»5#0x-  Hence  for  i  *  4.5  ve  find  2$  -  cosh  *  < 

Expression  (B.2)  becomes  for  t  i  4.5 

y  <  (l  +  y2  +  2y)  exp  (-  8 i)  . 

Integrating  this  expression  leads  to 

00  00 

[  Ex'  d*  <  |‘  E^  d«  <  |  (1  +  y2  +  2y)  e"36  . 

4.5  4.5 

For  the  range  y  <  15  in  vhich  E^  is  calculated  from  equation  (A.7),  ve  can  be 
assured  that 

B-l 


(B.?) 

-  36  =  -  8*. 

(B.3) 

(B.4) 
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j  Ex'  0*  <  1  x  ID"13  (  (B.5) 

4.5 

vhich  shows  integrating  to  $  ■  4.5  la  quite  sufficient  because  the  value  of  |  E^| 

In  this  range  of  y  Is  about  1  to  10.  Since  E^'  has  no  singularities  In  0  i.  ♦  <  4.5 
and  since  the  Integral  from  $  ■  4.5  to  infinity  is  finite  (and  very  small)  ve  can 
conclude  that  is  convergent. 

All  of  the  steps  after  expression  (B.l)  hold  for  Eg  as  well  as  E^.  Conse¬ 
quently,  the  inequality  (B. 5 )  holds  also  for  and  convergence  follows. 
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APPENDIX  C 

FORTRAN  IV  LISTING  OF  FTV  PROGRAM 


c  ####«  HTV  PROGRAM  «#*#* 

5URRCU  1  INK  PTV (TIMfcR2.T3»  T4*TG*RAO*PTS*OPTIUN,COSA»RHOW«CwaT. 

1  T  »  V  »  V  S  ) 

THIS  SUMPNOGHAM  CONTROLS  ThE  ITERATION  for  thf  PEAK 
translational  VELOCITY.  RTv.  IT  TS  The  only  SUBROUTINE  of  thf 
PTV  PROGRAM  which  is  CA|_LEn  From  the  MAIN  program. 

DIMENSION  OX(IOOO)  ♦OY(lnOO) *Is(P) 

01  MEN  S I  UN  G<6) 

DIMENSION  A (bO) .C (SO) 

COMMON  /UXYAJX.QY 
COMMON  /OIS/IS 


IF  (OPTION. GT  .0.)  GO  TO  1*) 

WRITE (fc.bao) 

WRITE (ft  *600 )  TIMERS*  T I*TA*Tb*RAO*PTS*OPTION*COSA.RHOW*CwAT 
WRITE (tb.bRO) 

74.P14S/  IS  A  ONUS  CoNvLRs I UN  FACTOR 
10  VC  =  2.* /4.21  4S7/RHCw/RA0 
N=PTS 
T=T4 

D T s ( T b -  I ) /FLOAT (N-l ) 

IE(T.LE .0. )  NaN-1 
IF(T.LE.O.)  T=DT/?. 

IS(1)=2 
IS(?)=I 

G ( ? ) *Cw  A  I /R A i ) 

6  (  3 )  s  T  I  M  E  R  ? 

G(S)=S0iRT  (TIMLR2-TJ) 

G  (ft ) =S0HT (TIMER?) 

C  INITIAL  SEARCH  FOR  MAXIMUM  Vt|OClTY 

DO  40  1=1 *N 
G(1)=T 
V  =  VC*E  V  (G) 

A  { I )  =  T 
C(i)=v 

VS=?.*COSA*V 

IF  (OPT  KJiM.Lt  .0.  )  WRlTF(ft»6lU)  T*V,VS 
T  =  T*0  T 


C-l 


n  n 
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4  (i  C«NT1*LF 

I  TLWAT  JO  J  (OP  PT  V 

Dt  TTh  ^  1  Nit  T  Hfc.  MAXIMUM  VfUJcITy  F  P(M  C  AKPAY 
CALI.  <Ma  (C*N*M,M1  ) 

A<?  =  A  (Ml) 

C2=c (w  1 ) 

A ( ] )sA (M) 

ciii*cm) 

A  (  ? )  r  A  ✓ 

C  (2)  s C.< 

OAsOl 

T=A  { 1  ) -1 •  M*U  A 

IF  (  T  » L E  •  0  •  }  f  =  l)A/5  « 

OT -\)k/?  • 

DO  45  J  =3 ♦ 1 U 
G(1  )  =  T 
VsvC«F  V <G) 

A  ( I )  s  7 
C  (  I )  =  V 

VS=2 • *COSA* V 

IF  ( OP  1  l O.M * L  t.O.)  «PI  It  (A«6)0)  T»V.V5 
T*T*m 
45  CONTI NLt 
(4=1  n 

IF  ( I  AMS  (M-m1  )  .LT..3)  UO  tD  Sb 
T  =  A  (?)  - 0 •  rt*D  A 
IF  (T.l  F  .0.  )  T=DA/5. 

UT*0A/3. 

DO  50  l  =  11.  lb 
D(1  >  =  T 
V  =  VC*F  v (0) 

A  ( T  )=T 
C  ( I )  =  v 

VS=?.«C05A*Y 

IF  (OPTION, Lt-.O.  )  <»WJ  TF  (M610)  T*V.VS 
T  =  T  ♦!)  r 
50  CONT  I '40t 
Ns  1  A 

S3  CONT I  Mjfc 

DO  75  J J= 1  *  b 

CALL  XV AX (C»N*M*M1  ) 

IF (JJ.LT. i)  00  TO  bp 

IF  (  AMS  (  (C  (m) -c  (Ml )  )  /C  (M)  )  .|  T.o. 001 )  00  TO  110 
IF(JJ.FO.b)  GO  TO  l?n 
h?  N=10 

T 1 sA { M ) 

T  2  =  A (Ml  ) 

Y1=C  <M) 

y/2=r  (."  I ) 

A (4)  =  I  I 
A  (  ]  0  )  =  J  2 
C  ( <4 )  =  Y  1 
C ( l n ) = v? 
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UT  =  APS(T1-1<2)/S, 

1 1  =  1 

DO  7  0  1  =  1 *  O 

T*T1 *01 *FLOAT ( < I-IO) /R*T 1 > 

IF  (T  .Lt.  0.0)  GO  TO  64 
G(1 )=T 

V  =  V C«F v (G) 

VS=£.*COSA*V 
GO  TO  6b 

c  WHEN  T  IS  LESS  THAW  ZERO  SET  TO  /FRO • 

64  T  s  0.0 

V  =  0.0 

60  IF  (OPTION  .LE  •  0.0)  wRjft  (6.610)  TtVtVS 
A  (  I  )  *  T 
C  (  T  )  =  V 
II=-1*II 
i'll  CONTINUE 
7S  CONTINUE 
110  V*C(M) 

T  =  A  ( M  ) 

VS=?.*COSA*C (m) 

IF (OF  I )ON,Lfc.n, )  WRlTF(6*b?0)  A(M)*C(M)tVS 
RETURN 
1/n  VsC(M> 

t=a (V) 

VS=?.*COSA*t (M) 

VS1=2.«C0SA*C(M1) 

WRITE (fc.6 TO)  T  tV*VS*A (Ml ) tC(M\ ) tOSl 
RETURN 
C 

c 

S‘<!)  Format  (  1H1  ,  10X»30HTRANSlATtONaL  VELOCITY  PROGRAM  ) 

SO.)  FORMAT  ( lH0,bX.4bHI  IERATtoW  FOp  PEAK  TRANSLATIONAL  VELOCITY  PTV  // 

1  1?X»RFT  1MK  (SlrC)  .HXtlbHvELoClTY  (FT/StC)  * 3X » 25HVERT I CAL  VELOCITy(F 
«2T/SFC)  /^yx«  I6HTAWGE1  SllHMFRGf  D*7x*17mTARGET  AT  SURFACE  ) 
r>00  FORMAT  (  lH0,bX.<i3Hlr)PUI  TO  SDUPODTTNE  PTV  //  lAX* 

1  45HT IMFH2, T J,T4»Tb»KAD,HTs»Upri0N*CUSA*RH0wtCwAT  //lPSEU.S/ 

2  IPSE  1<»  .5  ) 

610  FORMAT ( 1P3E22.6) 

620  FORMAT ( 1HO»6X*20H****«***************  , 9X » JHPT V « 1 9X * 3HPT V// 

1  1P3EP2.6) 

tHu  FORMAT ( 1H0.42H***  WARNING  ITERATION  DID  NOT  CONVERGE  «**  ,5X* 

1  3SHMA*IMU.-|  AND  NEAREST  VALUE  ARE  GIVEN  // 

1  1PX.RRTIM1  (SFC)»AX*l6HvtL0ClTY(FT/SEc)  »3X *2SHVERT ICAL  VELOCITyIE 

2  T /SEC  )  /***»  16HTAHGFT  SijMrtF  RG)  0  »  7  X  .  1  7hT  ARGE  T  AT  SURFACE  / 

3  (1P3E22.6)) 

C 

EN() 
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FUNCTION  F V ( 6 ) 

THIS  SUBPROGRAM  SETS  UP  THE  INTEGRATION  FOR 
THE  TRANSLATIONAL  VELOClTT  V 

DIMENSION  (i  ( 6 ) 

EXTERNAL  FI 
DATA  N/1H/ 

NNsFLOAT(N)*G(l)*G(2)/H, 

NNsMAXO (NN*H) 

NNxM I NO (NN*N) 

Xafi  { ]  )  ~8 •  /G  ( ?. ) 

IF ( X  »GT #  G  { J ) )  GO  TO  43 
ZlaG (H) 

IF(X.«U.O.)  Zl*SORT  (G(3)-X) 

Ir (0(1) • GT  «G ( 3 ) )  GO  TO  40 
G(4) s-1 .0 

Z2aSQRT ( G ( 3 ) -G (  I ) ) 

INTEGRATION  FOR  T  ,Lt ,  T1MFH2 
KVa-FGl  ( Zl «  ?2«NN»F 1 *G) 

RETURN 

40  z?=n . 

Z3aS0RT (G(l ) -G (3) ) 

IF(G(3) • L  0  «  0  •  )  GO  TO  4h 
G(4)*-l .0 

Nla71/(2l*/3)*FLOAT  (NNUtttQ 

NNNaZ.V ( Z 1 *Z3 ) *FLO AT ( NU ) +2.0 

INTEGRATION  EOR  INTERVAL  which  includes  timer? 

VU-FG1  (Z1*Z2,N1*F1*G) 

G ( 4 ) *  1 .0 

V2aFGI (Z2*73*NNN»F1*G) 

FVaVl *V2 
RETURN 

4 3  ZZaSQRT  (X-G ( 3)  ) 

Z3=SGRT (G(l ) -6 ( 3 ) ) 

<*o  G  (4)  *1  «U 

integration  fur  t  larger  Than  timer?  rut  the 

INTERVAL  OOtS  NOT  INCLUDE  TIMER?. 

KVaFGI  ( 22 • 2 J  t NM'F 1 «G ) 

RETURN 

end 
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FUNCTION  F 1 <Z,G) 

C 

C  THIS  SUBPROGRAM  CALCULATES  THf  PRODUCT  I NC I OE  N  T  PRESSURE  * 

REDUCED  STEP  WAVE  ACCELf«ATIOn  BY  CALLING  THE  INTERPOLATION 
PROGRAMS  vtab  and  ptah. 

DIMENSION  OX (1000) .QY(loOO) *IS(2) 

DIMENSION  G<6)  »QGIX  (120)  ,UQY<120> 

COMMON  /tilXY/QX  *QY 
COMMON  /OIS/IS 

c 

C  REDUCED  STEP  «lAV£  ACCELFRAtIOm  OF  A  CyLINDER 

C 

DATA  (OlM  (  1  ) , 1*1,106)  /O. » .0125, .026*. 0375, . 060 , . 075 » . 1 00  , 

1  •  1 ?5 • • lbO  « •  1 75* , 200  * . 2£5, , 25n , . 275, «  300  » • 325 . . 350 • • 37b, 

2  ,4n00* .425* .450* .A75..600, .5?b,  .SbO  t  .  S75  *  .  600  •  .625  *  .febO » 

3  .675, , 700.. 7?b,.750» .7  75, .HOo, ,025,.HSO* .875* ,900*.925,.950* 
A  .975* 1.0 0*1.05,1. 10*1 .1b,].2n*l.?5«1.30, 1.35. 1.40*1.45, 

b  1  .50*1  .55,  l«AO ,1.65*  1 .70*1 ,7s#  l  .0'),  1.85,  !  .90,1 .95*2 ,00, 

6  2. 05. c. 10.2. lb, 2. 20  ,2.?b,2.  Jo , 2 . 3b . 2 . 4 0 , 2 . 4b . 2 . 50 » 2 . 5b , 

7  2. 60. 2. 65. 2. 70. 2. 75, 2. 00*2. 8b. 2. 90. 3. 00* 3. 10. 3. 2. 3. 3. 3. 4, 

8  J,5*3.6*J.  7*3.0t3.9*4.O«4.2t4.4.4.6,4.P«b.O»6.?5*5.bO, 

9  5. 7b. 6. 00. 0.25*6,5.7.0*  7*6,8. 0  ✓ 

DATA  (LOY  (I)  *  1*1*60)  /  0.0,  .]  90 1  93 .  .275935  *.  33269A  ,.  37fl  1 80  • 

1  .4APHJ6, .SO2109..S44OOO* .b77l42*.60Al U  * .62bS«9, ,642701* 

2  .65614  3,  .466457,  .6  74 'i 7g,  .679  365 *  .6826 1 2  ,.  604  0 70  ,. 603955, 

J  .6H24S2, .679721 ,.675904, .671 127, .665499, .6591 2C.6520 70, 

4  .644453, .6303 15,.  627730* .610  755, ,609444* .599044, .509999, 

5  .579949*. 5 69730, • 559 3 ?4*,64rtQl3*. 5 30 372, .527777*. 517151, 

6  .506515* .495807,. 405204  * ,464215  * .4* 34 1 7 *. 422977 *. 402968* 

7  .  JP3447,.364460,.34t>ri4?,  .  320218,  .31 1008,  .294424,  .270471  , 

0  .2631b2,.24«465,.2344O4,.22O96O,.2O0l24,.195001  / 

DATA  (UUY( I) *1*61,106)  ✓  . 1 04219, . 1 73i22, . 162S73, . 152555, 

1  . 143051 , . 134041 ,. 125500, . 1 1 7435, .109flOl , . 1 02690 , . 095702 , 
c  .009361 ,.08 JlO0,.u 7 760flt .072242,. 06 7 196, .062453, .057999, 

3  .053818, .04  9097,. 0  4  6221  * • 03965b * . 033725 , . 0206 37 , . 024209 , 

4  .020 360,. Ol7044,,0 14  1  77*. 01171 2, .009699, .00  7795,.0 06260, 

5  .  0  03863*. 002 172,.  00 1 009, .000230,-0.00 0267,-0. 000619, 

6  -0,000774,-0. 000004,-0. 00 0  767 *-0.0  0  0696 *-0.00 0606, 

7  -0.000430,-0. 000297, -0,U0020b  / 

C 

TF (6(4) .GT  .0.  )  GO  TO  ?n 
X*G(3)-Z*Z 
GO  TC  30 
20  X-G (3 ) +  Z#Z 
30  XD* (G ( 1 ) -X ) *G (2) 

IF (7.GT • G ( 6 ) )  GO  TO  35 
PsPTAH  (X,(JX,UY*IS(2)  ) 

GO  TC  AO 

lb  P  =  VtAH<X»0X»0Y»IS{2)  ) 

4  0  F1=Z«P«VTAH  (XO,QQX,UUY ,TS( 1  )  ) 

RETURN 

end 
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SUBROUTINE  XMAX (B.N.M.Ml ) 


THIS  SUBPROGRAM  DETERMINES  THF  LOCATIONS  OF  THE  TWO  LARGEST 
ABSOLUTE  VALUES  OF  MEMBERS  OF  THE  B  ARRAY. 


DIMENSION  H(S0) 

X*AHS (B  (1)  ) 

Ms  1 

DO  10  1=2*N 

IF ( ABS ( B ( I ) ) .LT.X)  GO  To  lO 
Ms  I 

XsARS (H (M) ) 

10  CONTINUE 

Ml.l 

IF(M.EQ.l)  MIs2 
XrARS (B (Ml ) ) 

DO  ? 0  1=2. N 

IF  CABS  <B  ( I )  )  .LT.X)  GO  TO  20 
IF (I.EG.M)  GO  TO  20 
Ml  *  I 

X=A0S(b (Ml) ) 

20  CONTINUE 
RETURN 
END 
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FUNCTION  VTAB(X*Y,Z.K) 

THIS  SUBPROGRAM  PERFORMS  A  SECOND  ORDER  LAGRANGIAN  INTERPOLATION 

THE  INDEPENDENT  VARIARLF  IS  STORED  IN  THE  Y  ARRAY  IN  INCREASING 
ORDER.  THE  DEPENDENT  VARIABLE  IS  STORED  IN  THE  Z  ARRAY. 

X  IS  THE  POINT  AT  WHICH  THE  FUNCTION  IS  TO  Bt  EVALUATED. 

X  IS  Tht  NUMBER  OF  THE  ELEMENT  In  THE  Y  ARRAY  WHICH  Is  FIRST 
COMPARED  WITH  X. 


DIMENSION  Y  (1000) »Z (1000) 

IF (X.LE.0.)  GO  TO  SO 
UO  10  1 *K  « 1 00  0 
J*  I 

IF { Y ( I ) .GT.X)  GO  TO  20 
10  CONTINUE 
20  J*MAXO (3» J-l ) 

UO  30  1=1 f 1000 
IF(Y(J) .LT.X)  GO  TO  40 
J=J-1 

IF ( J.L I .3)  GO  TO  40 
30  CONTINUE 
40  K=J.l 

IF (? (J) .EO.Z (K) )  GO  TO  *0 
L*J-1 

A* ( X- Y ( K ) )/(Y(J)-Y(L) ) 

C=(X-Y (L) )/(Y(K)-Y (J) ) 

IF  c (A.LT.-S.O) .OR, (C.GT.S.O) )  GO  TO  60 
B«(X-Y(J> )/(Y(K)-Y(L> ) 
VTAR«C«(B*/(K)-A*Z(J))*A*B*Z(| ) 

RETURN 
SO  VTAH«0. 

RETURN 

60  VTARsZ ( J) ♦ (X-Y ( J) )*(Z(K)-Z(J))/(Y(K)-Y(J)) 
RETURN 
END 
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FUNCTION  PT  AH ( X « Y . 7 .K  ) 


THIS  SUBPROGRAM  PERFORMS  A  SECOND  ORDER  LAGRANGIAN  INTERPOLATION 
WITH  PROVISIONS  FOR  HANDLING  A  SINGULARITY. 

FUNCTION  ARGUMENTS  ARE  THE  SAME  As  IN  VTAB  . 


DIMENSION  Y (1000) * z  < 1 000 ) 

IFlX.LE.O.)  GO  TO  bO 
DO  10  IsK.1000 
J=I 

IF (Y ( I ) .GT.X)  GO  TO  20 
10  CONTINUE 
c'O  J=MAX0 (3* J-l ) 

DO  30  1*1.1000 

IF (Y (J) • L T • X )  GO  TO  40 

J=J-1 

IF { J.LT .3)  GO  TO  40 
TO  CONTINUE 
40  J=J*1 
JJsJ 

THE  FOLLOWING  THREE  STATEMENTS  PROVIDE  FOR  EXTRAPOLATION 
AROUND  A  SINGULARITY. 

IF (APS (7 ( J) ) .GT.1.0E2OJ jJ*J-2 
IF  { APS !/ < J-l)  ) .GT.1.0E20) JjaJ.l 

IF { tJJ.Eu.J) .AND. (ARS(Z(J-2> ) ,LT .  1 . 0E?0 ) )  JJ*J-1 

J=JJ 

Ksj.l 

IF (7 ( J) .E0.2 (K) )  GO  TO  60 
L  =  J-l 

A=(X-Y(K) ) / ( Y ( J) -Y ( L ) ) 

C* ( X- Y (L ) )/(Y(K)-Y(J) ) 

IF ( (A.LT.-b.O) .OR. (C.GT.b.O) )  GO  TO  60 

Hs (X-Y ( J)  ) / C  Y ( K ) -Y ( L ) ) 

PTAHsC* (B«Z (K) -A*7 (J) ) ♦  A  *B#2 (|  ) 

RETURN 
SO  PTAHsO. 

RETURN 

M>  PTAG*7 ( J) ♦ (X-Y  l  J) )*{/ (K)-7 { JU /( Y (K)-v ( J) ) 

RETURN 
t  NO 
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FUNCTION  F  6 1 (A*ti*K*F,P) 


this  subprogram  integrates  thf  function  f  betwfen 

A  ANO  8  USING  A  FOUR-POIhT  GAijSSIAN  QUADRATURE  IN 
K  SUBINTtRVALS. 


JO 

20 


DIMENSION  V<4) fW(2) ,SUM(4) eP(j  ) 

DATA. V/  -*«61 13631 1594053f«. 339981 043584856* 
1  .319981043584856,. 861136311594053  / 

DATA  W/  .347854845137454*. 652i4515486?546  / 

SUM(1)=0.0 

SUM(2)*o.O 

SUM(3)s0.0 


SUM(4) =0.0 


Hs  (8- A ) /E Lt)A T  (K) 

H2aH/2. 

AA«A*H2 

DO  ?0  L  =  1.K 

DO  10  1*1*4 

X*H2*V ( I ) . A A 

SUM ( I ) *SUM ( I ) +F (X«P) 

A AsA A  +  H 

FGI.H2«(W(1)#(SUM(1)4SUm(4))+W(2)*(SUm(2)*SUM(3) 


)  ) 
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APPENDIX  D 

SAMPLE  PROGRAM  OUTPUT  FOR  p(t)  «  exp  (-125t) 


translational  velocity  program 
Input  tc  subroutine  ptv 

TIMER^,T3.T4*TS»RADtPTS»0PTI0N*C0SA,RH0w.CwAT 


0. 

0.  7, 

B4000E-0?  2,20000E»01 

OOOOOE+OO  0. 

B.66030E-01  !* 

OOOOOE^OO  *5  • 00000E+03 

iteration  foh  pfak 

translational  velocity 

PTV 

TImE(SEC) 

VELCC1TY(FT/SEC> 

VERTICAL  VELOClTY(FT/SEC) 

TARGET  SUBMERGED 

target  at  surface 

9.800000E-03 

6.78763BE-03 

1.17S660E-02 

2.940000E-02 

7  *  1624  7  0t-04 

I.2405B3E-03 

4.900000E-02 

6.096742E-05 

1.055992E-04 

6.H60000E-02 

S.26O61OE-06 

9,11 1 692E-06 

3.<J20O00t-03 

6.2HB036t-03 

1.069126E-02 

i .372U0OE-0? 

4.8437B7E-03 

B.3B9729E-03 

2.352000E-02 

1  .S14464E-03 

2.623143E-03 

3.332000E-02 

4,347603t-04 

7.530309E-04 

4.312000E-02 

1.27149SE-04 

2.202306E-04 

5.292Q00E-02 

3.7349O3E-05 

6.469076E-05 

6*272 000E-02 

1.097070E-05 

1.900192E-05 

7.2S2O00E-02 

3.22202OE-O6 

5.5821 17E-06 

S.0S6000c;-03 

7*1  7db32t-03 

1.243365E-02 

1 .450400E-02 

4 • 463 1 B6L-03 

7.730S05E-03 

6*272 000E-03 

7.547924E-03 

1.307346E-02 

1 .332600F-02 

5.03915lt-03 

8.72BH2E-03 

7.44flo00l>03 

7.532047E-03 

1.304596E-02 

1 .215200E-02 

S.63756IE-03 

9.764594E-03 

8.624000E-03 

7.246194E-03 

I.255084E-02 

1 .097600E-0? 

6  »  233 1 34E-03 

1.079616E-02 

5.331200E-03 

7.289991E-03 

1.262670E-02 

7.212600E-03 

7.56029IE-03 

1.309488E-02 

5.566400E-03 

7.3B1641E-03 

1.278S45E-02 

6.977600E-03 

7,57b902E-03 

1.312365E-02 

5.801600E-03 

7.454637E-03 

1.291188E-02 

6.742400E-03 

7.5B077BE-03 

1.313036E-02 

6.036B00E-03 

7.S096S0E-03 

I.300716E-02 

6.507200E-03 

7  .  S7 1 3ij()t-03 

1.3I1403E-02 

#*#**#«*#*** «»«»«« « *  ptv  pjv 


6,742400E-03  /.bBi>77HE-03  1.313036E-02 
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